Galerkin Method in the Gravitational Collapse: a Dynamical System Approach 
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We study the general dynamics of the spherically symmetric gravitational collapse of a massless 
scalar field. We apply the Galerkin projection method to transform a system of partial differential 
equations into a set of ordinary differential equations for modal coefficients, after a convenient 
truncation procedure, largely applied to problems of turbulence. In the present case, we have 
generated a finite dynamical system that reproduces the essential features of the dynamics of the 
gravitational collapse, even for a lower order of truncation. Each initial condition in the space 
of modal coefficients corresponds to a well definite spatial distribution of scalar field. Numerical 
experiments with the dynamical system show that depending on the strength of the scalar field 
packet, the formation of black- holes or the dispersion of the scalar field leaving behind flat spacetime 
are the two main outcomes. We also found numerical evidence that between both asymptotic states, 
there is a critical solution represented by a limit cycle in the modal space with period An w 3.55. 
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Gravitational collapse is one of the most fascinating problems studied in General Relativity. Remarkable 
works helped the understanding of the process of gravitational collapse jl), but there are still crucial open 
questions concerning, for instance, the Cosmic Censorship Conjecture, and more recently the existence of 
critical phenomena in gravitational collapse. We are going to focus our attention here to this last issue. 

The pioneering analytical results of Christodoulou|^] together with numerical workjj| about the collapse of 
massless scalar fields basically established that, depending on the strength of the scalar field packet, two types 
of outcomes are expected to take place: the scalar field disperses to infinity leaving behind flat spacetime after 
an initial stage of collapse, or the scalar field collapses to form a black hole. Later ChoptuikQ], abandoning 
the restriction to almost stationary initial data, performed definite numerical work establishing the idea of 
critical behavior in the gravitational collapse and of a threshold of black hole formation. Originally, his 
numerical study was realized in the collapse of spherically symmetric distributions of massless scalar fields, 
and later several other fields such as perfect fluid 0, Yang-Mills^], charged and complex scalar fields 0, 
massive scalar fields||, gravitational waves [§, etc[lCj, were taken into consideration. Critical phenomena 
emerge as a robust feature in gravitational collapse, and might be attributed to the nonlinear nature of the 
Einstein field equations. 

Here we present some interesting preliminary results concerning the application, for the first time, of 
the Galerkin projection method[]ll| in the gravitational collapse of a massless scalar field. The idea is to 
treat the general gravitational collapse of a spherically symmetric massless scalar field without imposing any 
additional symmetry like some kind of self-similarity. The Galerkin method is an approximative technique 
that allows to transform a finite set of partial differential equations into an infinite countable set of coupled 
ordinary differential equations associated to the infinite modes of Fourier-type in which the solutions are 
decomposed. Truncation of these equations to a certain order of modes may give, under certain conditions, 
a finite dimensional dynamical system which approximates fairly well the exact dynamics. The application 
of this method in problems of fluid mechanics has brought interesting new results that has shed light in 
the understanding of turbulence, where the most famous application of the Galerkin method is the Lorenz 
system J53 (resulting from a truncation to third order of the Galerkin projection of Navier-Stokes equation) 
that describes the chaotic Rayleigh-Benard convection. 

Basically the Galerkin projection method consists in expanding a set of the system variables in a series of 
orthogonal polynomials (chosen as basis of the projection space X), substitution in the dynamical equations 
and obtaining a countable infinite set of ordinary differential equations for the modal coefficients. In order 
to carry out the Galerkin projection, the space X must be an inner product space and the basis functions 
must satisfy the boundary conditions implicit in the definition of X such that any candidate to a solution 
automatically satisfies the boundary conditions. Truncation of the system by setting modal terms beyond a 
certain order equal to zero (Galerkin approximation) results in a finite dimensional system which presumably 
yields an adequate approximation to the infinite dimensional dynamics if the truncation is at sufficiently 
high order. To gain insight in the types of dynamics that are possible we will consider here a low order 
truncation; while our truncation is not of high enough order to model the real behavior of our system in 
a completely faithful way, the resulting solutions give an indication of the type of qualitative behavior of 
which the actual system is capable. 

We consider the general spherically symmetric line element given by 



ds 2 = -gg du 2 - 2g du dr + r 2 dfl 2 , (I) 

where g = g(u,r), g = g(u,r) and dfl 2 = dO 2 + sinO 2 dip 2 ; u is the usual retarded null coordinate and r a 
radial coordinate which measures the proper area of the two sphere dfl 2 . The massless scalar field 0(it, r) is 
the only source of curvature, such that the relevant field equations that describe the dynamics are 

\$r (2) 

9 (3) 
v _ 

-j9A,r- (4) 

Eqs. (|2j) and (|3|) follow directly from the Einstein field equations, and the dynamics is completely encom- 
passed by the Klein-Gordon equation (Q). Another important quantity to be introduced is the mass function, 
m(u, r), 




3 



2m(u, r) 
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This quantity is interpreted as the effective gravitational mass inside the 2-sphere of radius r, and agrees 
with the Bondi and ADM masses in the asymptotic flat spacetimes. 

As we have discussed briefly the boundary conditions arc of fundamental importance for the Galerkin 
method, since they dictate the convenient set of basis functions. We then assume that the scalar field must 
satisfy the following boundary conditions, 4>(u, 0) = and </>(u,oo) = 0. Regularity of the spacetime at the 
origin requires g{u, 0) = 1, while the choice of u as the proper time at the origin fixes g(u, 0) = 1. 

After introducing the basic equations together with the appropriate boundary conditions, we are ready 
to apply the Galerkin method. For the sake of convenience, we introduce an alternative radial coordinate £ 
related to r by £ = Inr, where this new radial coordinate varies from — oo to +oo corresponding, respectively, 
to r = and r = oo. Then, the chosen orthogonal set of basis functions for the projection space X is 

MO = e-Z 2 / 2 H k (0, (6) 
where are the Hermite polynomials. The orthogonality is defined by the inner product 

/ + OC 
MO^j(Z)dt = 2 k ^k\5 kj . (7) 
-oo 

The following decomposition for the scalar field £) is proposed 



N 

0(u,O=X>k(u)tfk(O, (8) 

k=0 

where <pk{u) are the modal (Fourier-type) coefficients and N indicates the order of truncation. We remark 
that the precision of the above decomposition, in the sense of approaching the actual solution, depends 
directly on N. The boundary conditions for <fi are automatically satisfied since ipk(—oo) = ^(+oo) = for 
any k. 

The decomposition of <p(u, £) will be considered as the fundamental piece of our strategy of applying the 
Galerkin method to the system of Eqs. (0)-(|3). Before going through the Klein-Gordon equation, Eq. (Q), 
we need to express the metric functions g(u, £) and g(u, £) conveniently in terms of the basis functions. Thus, 
considering Eq. (En with the new spatial variable £, and introducing the decomposition of the scalar field 
given by Eq. dq), we have 



(9) 



The integration can be performed without difficulty for any order of the truncation; also from Eq. (^|) we 
have that g(u,t;) satisfies the boundary condition g — 1 at the origin £ = — oo, for all u. 

Now, the next step is to consider Eq. (^) that relates the metric functions g(u, £) and g(u, £). This 
equation can be written as 



9(u, = e-e J ^ a?g(u, O # = - |f + |§ - e"< / ^ || <, (10) 

where the term on the right-hand side arises from successive integration by parts. At this point, g(u,^) can 
be approximated in the following way 
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fe=i ^ 

Basically this approximation will be used in order to compute g from g, namely, as an integral of Eq. (flo|). 

The last step of the Galerkin method is to substitute Eqs. (||), (Q) and ( pi] ) into the Klein-Gordon 
equation that dictates the dynamics of the scalar field. The dynamical system for the modal coefficients 
4>k( u ) is obtained after projecting the Klein-Gordon equation into the nth mode ?/>n(£), n = 0, 1, 2, .., N, or 



fc=0 \ \ -s / / fe= Q 



^^(f^^))' » = 0,1,2..,JV (12) 



2 

fc=0 



where dot stands for derivative with respect to u. The result is a set of equations of the type 



Mu) =F k (<j> j ), fc = 0,l,2..JV (13) 

where the number of equations is defined by the order N of truncation. The expressions on the right-hand- 
side of Eqs. (P) are much involved, even for low N, and will not be given here. Once the modal coefficients 
are known either analytically or numerically, the overall spatio-temporal behavior of the scalar field and the 
metric functions are determined. 

Before the discussion of the qualitative aspects of the dynamical system, it will be very useful to 
establish a convenient way to choose the initial conditions. In our approach, a given set of initial 
conditions (0o(O), </>i(0), .., 0jv(O)) corresponds to the initial spatial distribution of scalar field given by 
0(0, £) = J2k=o 0fc(O)V'fc(£)- Thus, let us consider a Gaussian initial data 



0(0,0 =Ae 2 «e-^- ' 2 ) 2 , (14) 

where A is the amplitude of the distribution. The above expression can be decomposed with respect to the 
basis functions ifik{Q, determining the set of initial conditions, since 



such that the information about the initial strength of the scalar field is contained in the initial set 
(0o (0), 0i (0), ..,0at(O)), which depends directly on the amplitude A of the initial distribution of the scalar 
field. As illustration, for the value A = 2 of the amplitude of the initial distribution we obtain the initial con- 
ditions of the modal coefficients for the projections up to order N = 3, 0o(O) = 0.6782306563805649, 0i(O) = 
-0.05284561865739180, 2 (O) = -0.09123038744236496, 3 (O) = -0.0007815958656455472. We note that 
the higher the order N of truncation, the more the corresponding modal coefficients decrease; for instance, if 
we consider N = 9, 09(0) ~ 10~ 6 . A comparison of the original function ( |l4| ) with the initial condition func- 
tion constructed from the partial expansion of order N = 3 with the modal coefficients given above, shows 
clearly that the relevant contributions actually come from the lower order modes so that, in this aspect, our 
truncation to order N = 3 appears sufficient to grasp the essential features of initial scalar function wave 
packet. 

Now, we are ready to present the basic qualitative aspects of the dynamics in the space of variables <pk (u) ■ 
The origin is the unique critical point identified as the Minkowski or flat spacetime, because 0(«, £) = 
and g(u,£) = g(u,£) = 1. The linear analysis about the origin reveals that it is an attractor for any order 
N of the truncation. The eigenvalues associated with the linearized system about the origin are in general 
complex, conjugated in pairs, with negative real parts. From the physical point of view, this characteristic 
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is reasonable due to the fact that small modal coefficients fa(0) are equivalent to a small amplitude A of 
the initial data. As a matter of fact, the scalar field has not enough strength to maintain the collapse and it 
eventually disperses leaving behind flat spacetime. This behavior extends to a nonlinear neighborhood of the 
origin, generated from initial modal coefficients ^fe(O) associated with a finite (but not sufficiently large) value 
of the amplitude A. This behavior is illustrated in Figs. 1 for A — 2.0798, and is denoted subcritical. Here 
and in the remaining of the paper we are considering the Galerkin decomposition with TV = 3, meaning the 
presence of four modal coefficients. Fig. 1(a) shows the evolution of the modal coefficient fa(u); it oscillates 
for a certain range of u and then decays to zero. Other modes have a similar behavior. With the present 
initial distribution all modes tend to zero, namely, to the origin of modal space that represents flat spacetime 
(indeed the origin is an attractor for any order of truncation). Figs. l(b,c) depict, respectively, the scalar 
field 4> given by (|j) for N = 3 and the mass function m(u,£) defined from (||), as functions of £ for several 
values of increasing u. Initially <p and m may increase but as expected, for large u, both tend to zero, showing 
that the scalar field has not enough strength to collapse resulting in a final configuration corresponding to 
flat spacetime. The denomination subcritical given above is now clear, because of its analogy with similar 
behavior in the literature, in spite of our truncation at N = 3. Our computational resources allowed us to 
go to a truncation of order N — 5, with analogous results, suggesting that even a lower order truncation 
may grasp the skeleton of the full dynamics. 

A second possible configuration, denoted supercritical, corresponds to initial conditions connected to a 
distribution of scalar field with enough strength to hold the collapse until the formation of a black hole. The 
scalar field distribution in this case has amplitude A large than a critical value A* to be discussed later. 
The final outcome is the infinity region of the modal space characterized by X^fcLo ^fe = 00 ■ ^ n approaching 
this region the scalar field increases as well the metric functions g and g. Figs. 2 illustrate this behavior. 
Fig. 2(a) shows the modal coefficient <f>o(u) for the supercritical solution with A — 2.0799. After oscillating 
this mode increases in a very fast way. All remaining modes fa(u), fa(u) and fa(u) present an identical 
behavior. It can happen that at some point g/g — > 0, signalizing the formation of an apparent horizon. At 
this point the mass function is given by m a h — e^ ah /2 = r a h/2, with the subscript ah denoting the apparent 
horizon. Fig. 2(b) shows the expression 2m /r plotted as a function of £ for several values of increasing u. 
For u a h ~ 28.34, 2m/r w 1 at ^ a h ~ —0.54. This indicates the formation of a black hole. We remark that 
these numerical values are an approximation of those of the full dynamics due to the N = 3 truncation. 

So far the results are in agreement with previous numerical works ||. Following the terminology intro- 
duced by Choptuik, we called subcritical the first class of solutions, whereas the second class is known by 
supercritical. Between the two classes of solutions, there is the critical solution. The critical solution has an 
interesting symmetry known as discrete self-similarity, which is manifested by the periodicity of the metric 
and scalar field with respect to spatio-temporal scales. Gundlach and Koike et al|| [l(], engendered a 
very interesting picture in phase space for which the critical solution is represented by a limit cycle. In our 
dynamical system approach derived from the Galerkin method, we obtained evidence of a limit cycle in the 
space of modal coefficients, which corresponds to the critical solution constituting a limiting configuration 
between the subcritical and supercritical classes, according to the results of Choptuik. We proceed now to 
discuss our result. Consider A* the amplitude of the scalar field distribution associated with the critical 
solution. We have found that as A — > A* the modal coefficients oscillate for some time before escaping to 
infinity (black hole formation) or tending to the origin (flat spacetime). Numerically we can successively 
adjust A such that we approach A* from the region of initial data corresponding to black hole final state, or 
from the region corresponding to a flat spacetime final state. The more precise A is adjusted to approach 
the actual critical value A*, the longer is the time the modal coefficients oscillate periodically. Fig. 3(a) 
illustrates this behavior. The modal coefficients correspond to the choice A — 2.079856788846789, and it 
can be noted that after u rs 50 the oscillations are performed steadily with a definite period Au w 3.55. 
The Figure shows modes up to N = 3 but we have checked that all other modes, up to N — 5, exhibit a 
similar behavior also with Au « 3.55, indicating that we are approaching a periodic orbit in the complete 
modal space, or a limit cycle, which is a representation in the modal space of the critical solution. Fig. 
4 shows the projection of the previous solution in the 3-dim subspace spanned by the modes (fa, fa, fa). 
It is a remarkable fact that the presence of the limit cycle in the modal dynamical space manifests itself 
in any order of the truncation. This gives a strong indication that our dynamical system approach to the 
gravitational collapse, via the Galerkin projection method, exhibits the basic skeleton of the dynamics, as 
that which was found in numerical simulation of Einstein-Klein-Gordon system. We are presently making a 
careful evaluation of the order of the error between a truncation of order N and orders N + 1, N + 2, and 
so on. The indication is that the error tends to decrease fast as N increases, suggesting a faithful approach 
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of the dynamics of the system. 

We must finally remark the potentiality of the Galerkin method applied in the treatment of the dynamics 
of the gravitational collapse, as a reduction of the dynamics to a dynamical system problem, which may avoid 
extremely elaborated numerical techniques and, in its simplicity, guarantees that the resulting dynamics is 
a skeleton of the underlying physics and not a consequence of numerical artifacts. We are also presently 
applying the same method to the case of massive scalar field. 
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FIG. 1: (a) Evolution of the modal coefficient <f>o(u) for the subcritical configuration, which corresponds to the a scalar 
field distribution with initial amplitude A = 2.0798. The remaining modes have a similar behavior, (b) Behavior of 
the scalar field as a function of £ for several values of increasing u. The amplitude of <f>(£,u) increases indicating the 
initial phase of collapse, but for large u the scalar field tends to zero corresponding to a final flat configuration, (c) 
Distributions of the mass-function m(£, it) for increasing u. 
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FIG. 2: (a) Evolution of the modal coefficient <j>o(u) for the subcritical evolution with A = 2.0799. After an initial 
phase of oscillations, this mode increases in a very fast way. Other modes <j>\, (j>2 and <j>3 exhibit identical behavior, 
(b) 2m/r plotted as a function of £ for several values of increasing u. For this supercritical configuration, 2m/r ~ 1 
indicates the formation of an apparent horizon at {£ a h,u a h) and consequently a black hole. 
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First modal coefficient Second modal coefficient 




FIG. 3: Evolution of the modal coefficients <po : <t>i, 02 and (j>3 corresponding to the critical solution obtained with 
the choice A —= 2.079856788846789. Note that after u « 50 the oscillations are performed steadily with period 
Am m 3.55. By observing the amplitudes of the periodic motion along each modal coefficient, it is clear that they are 
smaller for higher order modal coefficients; for instance, for the first modal coefficient the amplitude is approximately 
0.8, and of order of 10~ 2 for the last modal coefficient used in our approximation. Eventually, if more terms were 
considered in the decomposition (8), it would be expected smaller corresponding amplitudes of the periodic motion 
of such modal coefficients. 




FIG. 4: The limit cycle present in the space of modal space for our case of truncation N = 3 projected into the 3 
dim subspace spanned by the modes (0o,0i,02). 



